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O Abstract 

In this work we present new second order semi-discrete central schemes for systems 
of hyperbolic conservation laws on curvilinear grids. Our methods generalise the 
two-dimensional central- upwind schemes developed by Kurganov and Tadmor [1]. 
c/3 In these schemes we account for area and volume changes in the numerical flux 

functions due to the non-cartesian geometries. In case of vectorial conservation 
laws we introduce a general prescription of the geometrical source terms valid for 
various orthogonal curvilinear coordinate systems. The methods are applied to the 
two-dimensional Euler equations of inviscid gas dynamics with and without angular 
momentum transport. In the latter case we introduce a new test problem to examine 
the detailed conservation of specific angular momentum. 
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1 Introduction 



In the last decades various multidimensional numerical schemes for the so- 
lution of hyperbolic conservation laws have been developed. The continuous 
increase in computing power provides the opportunity to perform accurate 
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computations even in the case of three-dimensional problems. However for 
some physical systems this might not be the appropriate approach. In cases 
where the underlying physics exhibit some kind of symmetry more accuracy 
can be achieved by omitting the dependence on one dimension. In some cases 
these symmetries might lead to different conservation laws. For instance in 
rotationally symmetric inviscid fluids angular momentum is locally conserved. 

Among the variety of numerical methods the family of Godunov-type schemes 
are a very useful approach. Since the seminal work of Godunov [2] these 
schemes became an important tool in the numerical analysis of hyperbolic 
conservations laws. The original proposition of Godunov was quite simple: 
Approximate the initial condition by piecewise constant data. Then advance 
the solution in time by solving the intermittent local Riemann problems. Since 
those days the method was improved by introducing less time-consuming 
approximate Riemann-solvers. Furthermore the development of higher-order 
methods yields better convergence of the numerical schemes. In most parts 
of our report we will follow the work of Kurganov and Tadmor [3|T] . They 
proposed Riemann- solvers-free second order methods which avoid computa- 
tion of the eigensystem of the advection problem. We incorporate orthogonal 
curvilinear coordinate systems into their cartesian scheme and therefore allow 
for area and volume changes of the grid cells. 



The outline of our paper is as follows: In Section IT we briefly review the co- 
variant formulation of conservation laws and derive the basic two-dimensional 
integral equation for general orthogonal coordinates. This result is utilised 
in Section [2] to obtain second order semi-discrete non-oscillating numerical 
schemes resembling those described in [T] . In Section [3] we present the results 
of numerical computations obtained with the new scheme for the equations of 
gas dynamics. The simulations were carried out on curvilinear meshes and on a 
two-dimensional cartesian mesh for comparison. These tests examine the solu- 
tion of two-dimensional Riemann problems in polar and cylindrical symmetry. 
In addition to that we introduce a method for testing the detailed conservation 
of angular momentum for rotationally symmetric flows in Section 3.3| Finally 
a summary is given in Section |4} 



1.1 Conservation laws and orthogonal curvilinear coordinate systems 



The concept of conservation is fundamental to a variety of physical phenomena 
and leads to partial differential equations of almost the same kind in very 
different areas. In this work we will focus on systems of non-linear hyperbolic 
conservation laws of the form 

^ + V-T(u) = 0. (1.1) 



Here u is either a scalar or vector and T a vector or rank-2 tensor. V- denotes 
the covariant derivative with respect to some affine connection followed by 
a contraction over the last indices. In the following we will restrict ourselves 
to orthogonal curvilinear coordinates {£, r), 0} with a local orthonormal basis 
{et, effi ei} and metric scale factors {h^, h v , h<f,}. For the divergence of vector 
fields one obtains for all components with respect to local orthonormal frames 
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and in case of tensor fields 



V-T 



V-T 



V-T 



)g (| ( Wff) + | (V ,%) + A (h,h v T^ 



~ c ivi T H ~ c 4>v4> T 4>4> + c viv T iv + c v4>v T < 
~ c iM ii ~ c v4>v w "^ c 4>i4> i<t> ~^ c 4>v4> v 




1.3) 



In these equations new geometrical quantities arise: The square root of the 
metric determinant y^ = h^h v h<f, and the commutator coefficients Cyfe which 
depend on the metric scale factors and their derivatives. A more detailed 
derivation is given in the appendix. 



In this paper we will focus on two-dimensional conservation laws by assum- 
ing a coordinate symmetry with respect to 0. We may consider the 3D case 
in a follow-up paper. Hence we demand that all functions - geometrical scale 
factors as well as physical quantities - are independent of 0. Therefore the com- 
mutator coefficients with in their second index vanish and a two-dimensional 
conservation law is obtained for the scalar variable u and vector field v(u) 
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In the same way we derive a vectorial conservation law for the vector w and 



tensor field T(w) 
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The only differences between scalar and vectorial conservation laws are the 
geometrical source terms in case of the latter. Therefore we can combine both 
to a system of conservation laws. At this point it is convenient to define new 
spatial differential operators 

V, = j_^h v h„ V n =±=^hth+ (1.6) 

and rewrite the conservation law 

d t u + V^F(u)+V v G(u) = S(u). (1.7) 

In this compact notation u denotes the vector of conservative variables, F(u), 
G(u) and S(u) are the flux vectors and geometrical source terms, respectively. 



To allow for discontinuous solutions one integrates the differential equation (1.7) 
over the time interval [t n ,t n +i] and a spatial region D given by the cartesian 
product [£_ , £ + ] x [r/_ , r/ + ] . Hence we obtain an integral equation, the so called 
weak formulation 



tn+l 



(u(t n+ i)) D = (u(t n )) D - J (V^F + V r ,G - S) D dt (1.8) 

t n 

using the notation ( ) for volumdM averages as defined in 
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1 The integration with respect to <\> is suppressed throughout the whole paper, 
because all functions are considered independent of (j). Nevertheless we will stick to 
the term "volume" for integrals over two-dimensional domains. 



with spatial volume AV of region D 
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Equation (1.8) describes the time evolution of volume averaged conservative 
variables u in region D. So far it is impossible to evaluate the flux integrals 
without further knowledge of the function u(t, £, rj) on the surface of D and at 
time t within the interval [t n , £ n +i]- However in the next section we will derive 
a numerical scheme which provides an approximation to these integrals. 



2 Numerical scheme 



2.1 Semi-discrete scheme for generalised orthogonal coordinates 



The derivation of the numerical scheme follows the three steps of recon- 
struction, evolution and projection described in pp. For illustration consider 
the control volume selected by the cartesian product D 
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in curvilinear orthogonal coordinates shown in Figure 



There 



are two types of staggered grid cells drawn along the boundary: 



2.1 



Edge cells are defined by the set union of partial regions at the edge of two 



adjacent cells (light grey in Fig. 2.1), e. g. at the eastern boundary: 



D l+y = Dl 3 UDr +hj 



whereas the partial areas addressed by the index pair {i,j} are given by 
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Corner cells are formed by the partial regions of the four neighboring cells 



which meet in a cells corner (dark grey in Fig. 2.1), e. g. around the south- 
eastern corner: 
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With the help of the definition of all partial corner areas with respect to 
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it is possible to construct the staggered corner cells. 
Central region: Finally the remaining central region of the cell is defined by 
the complement 
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Fig. 2.1. Schematic view of the control volume 

For the derivation of the numerical scheme we assume that at time t n volume 
averages of the conservative variables «"• := (u(t n )) D . are available for each 
cell in the computational domain. Cell boundary data is obtained via piecewise 
linear reconstruction. Using cell mean values and proper approximations for 
the slopes the linear expansion yields 
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. It is essential to deplete arti- 



ficial oscillations caused by this linear reconstruction process in order to obtain 



stability for second order numerical schemes. Various methods are discussed in 
the literature to achieve stable non-oscillating second order schemes, see e. g. 
jlfofBl?] . The original scheme by Kurganov and Tadmor follows the MUSCL 
(Monotone Upstream- centred Schemes for Conservation Laws) approach first 
proposed by van Leer [8]. This method introduces non-linear functions - so 
called (slope) limiters - to damp spurious oscillations. The method may also 
apply to a curvilinear scheme, if one demands consistency with the averaging 
process 






{ulj)^ = (u(t n )) D .. = ul 3 (2.3) 

This equation should hold independently of the actual choice for the slopes. 
Monchmeyer and Miiller [9] showed that this property is essential to retain a 
conservative scheme in case of non-cartesian grids (see also [10J ) . With help 



of Eq. (2.3) one derives the corollary that the coordinate pair (£o>?7o) is deter- 



mined by the barycentre of each control volume 

£o = (0 DiiJ Vo = (rf) Di . ■ (2.4) 

Therefore cell mean values might be regarded as point values at the barycentre. 



Equation (2.2) together with an integral conservation law for all the staggered 



grid cells similar to (1.8) forms the building block to advance the solution in 
time. As an intermediate result one computes cell mean values at time t n+ i 
defined on the staggered grid. Finally the updated mean values are obtained 
by reconstructing the staggered data and projecting these functions onto the 
original cell area D, h j 

<r = ^r_j^ +1 (tn)AV = {^) Dii (2.5) 

This is the generalised formulation of equation (2.3) in pQ for orthogonal 
curvilinear coordinates. The function w n+1 {^ 1 rf) refers to the combination of 
piecewise linear reconstructions on the staggered grid. Depending on the con- 
trol volume in which they are defined this function may vary from cell to 
cell. Nevertheless one can subdivide the integration into parts carried out over 
regions determined by the intersection of D^ with the staggered cells (cf. 



Fig. 2.1) 
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Up to this point no information about the advection problem has entered our 
considerations. This changes if we fix the limits for integration by introduc- 
ing the minimal and maximal local speeds of propagation for discontinuities 



according to [TT], e. g. at the western and eastern cell boundaries 
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Here C denotes a curve in phase space connecting two adjacent states uf_i 
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and u- i • of neighboring cells via the Riemann fan and A m j n , A max the minimal 
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and maximal eigenvalue of the Jacobian (^f-)- Similar definitions apply to the 



\ du 

southern bf-_% and northern &*_•, i wave speeds with the exception that the 

« 2 « + 2 

Jacobian of F has to be replaced by (f^)- Although these expressions seem 
to be difficult to evaluate, for genuinely nonlinear or linearly degenerate waves 
it is sufficient to compute (cf. [llj) 
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This is strictly true only for cartesian coordinates. The reason for this lim- 
itation is that the underlying formalism is based on the characteristic de- 
composition of the quasi-linear conservation law. In cartesian coordinates the 
transformation to a quasi-linear form is straightforward. A simple calculation 
with help of the chain rule yields 

du (dF\du (dG\du _ 
dt \du ) dx \du J dy 



However the curvilinear advection operators (1.6) involve derivatives of geo- 



metrical scale factors and the quasi-linear conservation law becomes 



9tu + (jj^JVp + { ^ W = S (u). (2.9) 



du 



This is equivalent to Eq. (1.7) if and only if the flux functions F and G are 
homogeneous functions of the conservative variables u. For a homogeneous 
function 



holds (cf. [12] Chap. 16.2). Hence 
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and (1.7) may be rewritten in the quasi-linear form (2.9). Apart from this 
discrepancy there is another pitfall when using curvilinear coordinates. The 
extent of the staggered grid cells is computed via multiplication of the local 
propagation speeds with the time step At. However in curvilinear coordinates 
spatial distances are not equal to coordinate distances. Therefore proper dis- 
tances should be obtained via evaluation of path integrals in compliance with 
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The accuracy of the approximations is sufficient as long as the coordinate 
distance is small enough. In fact in the limit At — > the staggered grid cells 
will collapse to lines, so that the considerations will hold. Hence the limits of 
staggered zones along the edges with respect to cell D it j are given by 
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Whereas in the corners one computes 
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Here the propagation speeds are derived from the values of neighboring cells 
according to (see [TT] ) 
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With these definitions we may expand the integrals over staggered cells arising 
in (2.6), e. g. for the southern domain D\ • the integral of an arbitrary function 

f(€v) is S iven by 
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With the help of fl2.1ip this leads to 



J f(Z, v ) dV = btj_iAt J mvJhK d£ +0(At 2 ). (2.14) 

Therefore the volume of the southern region may be expressed by 



AV% =/ dV = b+ j _ h te]h(h4d£ 



+ At 



^-i 



(2.15) 



10 



Furthermore it is necessary to compute approximations for the flux integrals 



arising in (1.8), e. g. again for the southern domain this yields 
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Similar equations can be obtained for the other staggered domains along the 
edges. To proceed with the derivation of the numerical scheme we expand the 



first sum in (2.6) up to first order in At using (2.14, 2.15) 
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where w™ +1 denote the staggered reconstructions. They are defined in the 
same way as the non-staggered reconstructions, e. g. in the southern domain 
according to 
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with mean value wfj x in the staggered area D: ,-_i = £),",• -, UD-,-. We would 
like to emphasise that the cell barycentres (£oj?7o) hi equation 2.19 are not 



identical to those in equation |2.4 because they depend on the domain under 
consideration. As in the case of non-staggered reconstructions we demand 



consistency with the averaging process [see Eq. (2.3)], i. e. 
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In case of the corner regions (second sum in 2.6) one can avoid the detailed 



computations. A short calculation proves with help of Eq. (2.12) that all vol- 



ume elements are of order At 2 (cf. [TTJ), e. g. for the north-western domain 
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Therefore using piecewise lineare reconstructions similar to those in (2.19) one 
concludes that 
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Henceforth one proceeds with the simplification of Eq. (2.6). The expressions 



(2.1812.20) replace the first and second sum of (2.6) and the updated cell mean 



values become 
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Here we used the constraint (2.3) to substitute (w™ + ) D5 by the mean value 



w™/ 1 within the central region. The next step in the derivation of the update 



formula for cell mean values incorporates the integral form of the conservation 



law. With help of (1.8) one replaces the mean values on the staggered grid at 



time step t n+ \ with flux and source term integrals. The integrals with respect 
to time may then be approximated by the midpoint quadrature rule. Hence 
for the central region one obtains 
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To simplify integration over the irregular shaped domain D\- one can substi- 
tute this integral with that over the whole cell area Dy and subtract those 



integrals over the supplementary domains along the cell boundary (see 2.1 ) 
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Again it is safe to omit the integrals over the corner areas (see 2.20). Since 
some of the flux integrals over edge domains are of order At (cf. Eq. 2.16.2.17), 
multiplication with At again leads to terms of order At 2 . Thus dropping all 
terms of order At 2 the contribution due to the central region reduces to 
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This result is completely determined by reconstructed data inside Dij whereas 
for the staggered cells around the boundaries different reconstructions of ad- 
jacent cells have to be taken into account. The same argument regarding flux 
integrals over boundary areas as mentioned above applies in this case. Thus 
the cell mean values within the edge regions are determined by 
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Insertion into (2.21) yields 
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Utilising (2.14) and (2.15) it is straightforward to prove that the contribution 



of all terms inside the curly brackets is of order At. Therefore one can divide 
the whole equation by At and compute the limit At — > 0. For convenience we 
define the time-dependent numerical fluxes across the cell boundaries 
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at the eastern edge and 
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at the northern edge respectively. Using (2.14), (2.15), (2.17) and similar equa- 



tions for other boundaries one derives the limits of all terms within the curly 



brackets in (2.22) 
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The remaining integrals on the left hand side of these equations cancel out 



with the cell mean values of flux derivatives in (2.22). Finally we reach the 



semi-discrete update formula for the volume averaged conservative variables 
within region D^j 
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with numerical flux functions given by (2.23) and (2.24). Cell mean values of 



source terms S and volume elements AVy should be obtained by integration 
over the domain Dij according to (1.9) and (1.10). 



Compared to the semi-discrete equation for time-evolution derived in pQ our 
result offers a higher degree of generality in two ways: 

(1) The scheme applies to orthogonal coordinate systems and therefore ac- 
counts for area and volume changes of grid zones. 

(2) It utilises integral representations for the numerical fluxes instead of ap- 
proximations with quadrature rules. 
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The latter point allows us to discretise the numerical fluxes and source terms 
in different ways starting from the same integral formulas. 



2.2 Numerical approximations of flux and source term integrals 



The accuracy of the numerical scheme described in the previous section is 
limited by the order of the polynomials in the reconstruction process. Hence 
using quadrature rules of higher order to approximate the flux and source 
term integrals do not improve the overall accuracy of the scheme. Whereas 
quadrature rules of lower order would compromise the second order accuracy 
of the scheme and lead to a larger numerical dissipation. Therefore we will 
focus on second order quadrature schemes like the midpoint and trapezoidal 
rule. 

In the former case we define the reconstructed values at the cell interfaces 
according to 



u 



hJ ~- Uij (&_ 1,1] j) Ulj = U id (& , T]j_ i ) 

u lj = &i,j (&+ 1 , Vj) K,j = Ui,j (& , Vj+ 1 ) 

and the numerical fluxes are given by 



(2.26) 



-<kJ a »U<i-<-v)\ (2 - 27) 



3 + 1 J 3+1,3 < 



with area elements being 



A^i+lJ — hrjhq 



- b l + y b inM3-<3+j} ( 2 - 28 ) 



A v , AA i>j+ i = hzht Af . (2.29) 



Z i+ i,Vi " 



The same rule applied to volume elements and source terms yields 

AV™< = hzhnhj, c m A£Ar) (2.30) 



liWi 



( S )dI = -XTT^ SiUi^htW^AZATi = Sfakm). (2.31) 



3 A K 
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In case of the trapezoidal rule we have to define the corner values as 



u 



Z = U iMi-\^3-\) 



M M 



UiMi+hVj-k) 



u 



ij = Ui J (&- 1 ' %+ h ) Kj = Ui J (&+ £ ' %+ 1 ) 



(2.32) 



and obtain 

-rtr 

vT-,1 — 



2 fe 



— 0- , 1 ■ 



a4 + , j+ , Ui.^K-) - <t f KTi j} 



■!,(' 



- a i+ij a i+lj^i -W i+lj 



<y-iK^i f ^) ( 2 - 33 ) 



^/(C 



se sw 

a i+l, j a i+hj\ U i,j ~ U i+l,j 



OZ-i 



3+ 



1 

' 2(6+ -6- , 



A 4^ +i fo i + iG«)-^ +i G(^ +1 ) 



6 &^ w<5 - «5 + + A ^t*,i +4 & i+^«i 



(2.34) 



-^ + ^(< + i)-ft + |^ + i(<i-< + i)^ 

for the fluxes. Here the area elements are given by 

AAf +i i = h v h^ At], AAf +i x = h v h 



Ca.Ifl 



»+o J+o 






A4^ 



*+h3+i 



h^h q 



^'"s+h 



a& a4_,. +4 =^ 






A77 (2.35) 
Af. (2.36) 



The volume elements yield 



rtr _ 


\ h h h \ h h h 


»J 


\ l 2 J 2 l+ 2 J 2 




+hch„h ( j ) + hfh„hs 



and source terms, respectively 

(S) 



mr 

'a, 



- S h^hnhtj 



4 AK- 



*,3 



t_l »»?,•_ 1 



+ 5 h^hrjhq 



-S h^hrjhq 






o hch^rLA, 



t i+ ^j-i 



^J'VJ, 



(2.37) 



(2.38) 



In the cartesian limit all scale factors are identical to one and these formulas 
reduce to those derived in pQ. 
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If the conservation law under consideration has a vectorial form, inertial forces 



appear in the curvilinear description (cf. Eq. 1.5). To account for these geo- 
metrical source terms we introduced the concept of commutator coefficients in 
Section [TTTj These non-linear functions depend on the scale factors h^, h v , h<f, 
and their derivatives with respect to the coordinates. Numerical approxima- 



tions of the commutator coefficients may in principle involve Eqn. (2.31) and 



(2.38). The remaining question is, how to approximate the derivatives of the 
scale factors. A comparison of the truncation error of the source term integrals 
with the flux differences shows that in case of the midpoint rule central dif- 
ferences perform better whereas for the trapezoidal rule one-sided differences 
lead to better results. Therefore the commutator coefficients for the midpoint 
rule are given by 



Cr 



Ar) 



c 'vh = ~ A^M&fy) W&+£>»fe) -hn{ti-$>Vj 









cZ 



(2.39) 






h3 



In case of the trapezoidal rule one requires four corner values for each of the 
four commutator coefficients, i.e. 



C ^ ^^H^HHilM^^HJ-M^-^H 



^tr,sw 

~ A AT, 

tr,se = 1 AT] 

C vifi 4 AV i>3 



';j V 

4 \ 



^tr,se = 1 A?7 

'<#} ~ AAV,, 



(2.40) 



and similar expressions for the remaining commutator coefficients cg-g, ci 
and cxaI- 



3 Numerical Experiments 



In the numerical examples we focus on the Euler equations for inviscid com- 
pressible and non-heat-conducting gas dynamics. This system of non-linear 
conservation laws may be written in generalised orthogonal coordinates ac- 



cording to (1.7). Thereby we assumed symmetry of the flow with respect to 



the spatial coordinate <fi. The vectors of conservative variables and fluxes in 
the £ and r\ directions are given by 



it 



Q 

QVf) 

Q v $ 

E 



gv 2 ~ + p 

QVgVr) 

(E + p) Vi 



G 



QVf, 

QVf,V^ 
QV?+p 

QVf) v 4> 
(E + p)vfi 



(3.1) 



If we furthermore assume, that the equation of state is determined by the ideal 
gas equation, the total energy E depends on the density g, the pressure p as 
well as the velocity components vt, vq and vi, respectively 



E 



~ e (yl + $ + >$) + 



p 



7 



(3.2) 



Here 7 denotes the constant ratio of specific heats. It is set to 7 = 1.4 in all 
test configurations. For the geometrical source terms one computes with the 



commutator coefficients (A7) 
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gv£ 


















~ Vf i c zf)i 
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C U4> 




v i c M 


+ QVf, 


~ v £ c viv 


+ Qv\ 


4>f)4> 


+ p 


~ V 4> C U4> 




~ v 4> c 4m4> 

























c v£v + c 4>£4> 



c^ + c 








(3.3) 



Besides these geometrical sources we do not account for additional body forces. 
Therefore all units can be removed from the system given above. In all tests 
we use this dimensionless prescription of the equations of gas dynamics. In 



Sections 3^ and 3^ we examine Riemann problems with the initial data of 
vi set to zero. In this case vi remains zero and the fourth component of the 



vectors in (3.1) and (3.3) vanishes. The system of Euler equations reduces to 
a pure two-dimensional problem without any flow in the direction of <fi. In 



Section |3.3| we study three-dimensional flows with rotational symmetry. The 
fourth equation of the above mentioned system then describes the conservation 
of the angular momentum. 



To advance the numerical solution in time according to (2.25) we used a third 



order Runge-Kutta scheme as described in [T31IT4"] . Cell boundary data is ob- 
tained via piecewise linear reconstruction. Unless stated otherwise we apply 
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the monotonized- central limiter proposed in [15] 

mc (A H ,,, A i+hj , A itj ) = minmod (OA^y, A itj ,6A i+ i d ) (3.4) 

with parameter 9 G [1,2]. The left-handed, right-handed and central differ- 
ences - e. g. in the ^-directions - are given by 

and the multivariable minmod function is evaluated according to 

{mm(xi, . . . ,x n ) ifXi>0Wi, 
max(xi, . . . ,x n ) iixi<0Wi, (3.5) 

otherwise. 

The amount of numerical diffusion is controlled by the parameter 9. Lower 
values are more diffusive while higher values result in an unstable scheme. In 
most of the test problems we set the parameter 9 = 1.3. This is less diffusive 
than the ordinary minmod limiter (which corresponds to 9 = 1.0) but retains 
stability in most cases. 



3.1 Two-dimensional Riemann problem on a polar mesh 



In this section we analyse the numerical solutions of two-dimensional Riemann 
problems on a polar mesh. The same tests were used by Kurganov and Tadmor 
PQ for comparison with the numerical results of Schulz-Rinne et al. [16] as well 
as Lax and Liu [17] . The computational domain is a square in the cartesian 
case and a circle for the polar mesh with the origin in the centre in both cases. 
At time t = the simulation is initialised with piecewise constant data in the 
quadrants defined by the x and y axis. 

Kurganov and Tadmor studied the numerical solutions of 19 different config- 
urations. Since our numerical schemes differ only from those described in [1] 
by means of the geometrical factors, we will focus on the discrepancies due 
to the polar mesh. Differences between computations on cartesian and polar 
grids would appear in each of the 19 test cases. Therefore we show only our 
results obtained for test configuration number 18 in [TJ with the initial data 



Vi = i 


£2 = 2 


Pi = 1 


Ql = 1 


u 2 = 


v 2 = -0.3 


m = o 


Ul = l 


p 3 = 0.4 


£ 3 = 1-0625 


p A = 0.4 


Qi = 0.5197 


% = 


v 3 = 0.2145 


M 4 = 


v 4 = 0.2741. 
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Fig. 3.1. Density at time t = 0.2 computed on a cartesian mesh (left) and polar 
mesh (right) using midpoint rule for flux calculations. 30 equally spaced contour 
levels between 0.525 and 2.025. 

u and v denote the velocity in the x and y direction, respectively. The north- 
eastern quadrant has index 1, the others are labeled counterclockwise in as- 
cending order. These initial conditions generate a shock-wave between quad- 
rants 2 and 3, a rarefaction- wave between quadrants 1 and 4 and contact 
discontinuities between quadrants 1 and 2 as well as 3 and 4. In the centre of 
the computational domain the 4 different solutions join leading to a complex 
flow. 



The system of equations under consideration is given by (3.1), (3.2) and (3.3) 
with vi set to zero. For cartesian coordinates all scale factors are unity and 
the commutator coefficients and thus all geometrical source terms vanish. In 
case of the polar grid we identify £ = r, rj — <p, <fi — z, respectively and assume 
slab symmetry, i. e. all derivatives with respect to z are zero. The scale factors 



are h r = 1, ha 



and h z = 1. The remaining non- vanishing commutator 



coefficient is therefore 



viv 



^ipTLp 



1 dh v 
hrh^ dr 



The computational domain covers a region of 1 x 1 in non-dimensional units 
with a resolution of 400 x 400 cells on the cartesian mesh. In case of the 
polar coordinates the extent of the circular domain is so large, that the unit 
square fits exactly into it. Thus the radius of the computational domain is 
\/2/2 with a resolution of 282 cells. The angular resolution is 360 for the full 
27r circular domain. The parameter for the limiter was set to 9 = 1.3 in all 
examples and the Courant number is 0.4. A stability criterion (CFL condition) 
similar to that proposed by Courant, Friedrichs and Lewy [T5] for an explicit 
central-upwind scheme is given in [3J. 
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Fig. 3.2. Density at time t = 0.2 computed on a cartesian mesh (left) and polar 
mesh (right) using trapezoidal rule for flux calculations. 30 equally spaced contour 
levels between 0.525 and 2.025. 



Figures 3^ and 3^ show the density contours for this two-dimensional Rie- 
mann problem. The numerical fluxes were calculated utilising the midpoint 



rule (Fig. 3.1) and the trapezoidal rule (Fig. 3.2). There is almost no visible 
difference in the dependency of the numerical fluxes. A quantitative analysis 
shows relative deviations of the density distributions of less than 10~ 3 except 
for the innermost area. Within a region of roughly 0.1 around the centre one 
observes a discrepancy between the solutions obtained with the two different 
flux functions of a few percent. The results given in [T] for cartesian geometry 
are similar to ours besides the small 'bump' in the north-western quadrant. 
In our results this 'bump' shows up only on small scales at a density level of 
£ = 2.0125. 

In the simulations on the polar mesh we achieved results which are almost 
identical to the cartesian case, at least in the central region. At a distance of 
0.2 from the centre the resolution of the shocks and contacts is less sharp on 
the polar mesh and they tend to fan out towards the boundaries. The same 
applies to the rarefactions in the north-eastern quadrant. The main reason for 
this may be the decrease of the angular resolution with increasing distance to 
the centre. 



3.2 Spherical Riemann problem between walls on a cylindrical mesh 



Langseth and LeVeque proposed a spherical Riemann problem in [19]. The 
flow under investigation is rotationally symmetric and may be examined using 
cylindrical coordinates. Therefore we identify £ = z,r] = r,(f) = (p and assume 
symmetry with respect to the polar angle if. The geometrical scale factors 
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Fig. 3.3. Pressure in the r-z plane at t = 0.7 computed with different numerical 
fluxes using the midpoint rule (left) and the trapezoidal rule (right). 32 equally 
spaced contour levels between 0.73 and 1.48. 




Fig. 3.4. Horizontal profile of the pressure at z = 0.4 and time t = 0.7 computed 
with different numerical fluxes using the midpoint rule (left) and the trapezoidal 
rule (right). 

are given by h z = 1, h r = 1 and h v = r. In this case the only non- vanishing 



commutator coefficient is 



CZf.2 — C, 



iprip 



1 dh v 

hphr dr 



The initial condition is a homogeneous density distribution go = 1 with van- 
ishing velocities. Inside a sphere with a radius of r = 0.2 centred at z = 0.4 
the pressure is set to 5 and outside to 1. The fourth equation can be removed 



from the system given by (3.1), (3.2) and (3.3), because initially there is no 



angular motion. Hence v v remains zero. The computational domain covers a 
region of [0, 1.0] x [0, 1.5] in the r-z-plane with a resolution of 400 x 600 cells. 
The parameter for the monotonized- central limiter is set to 9 = 1.3 and the 
Courant number is 0.4 in all computations. 



In Fig. 3.3 the pressure at time t = 0.7 in the r-z-plane is shown for the two 
different numerical fluxes. Both methods seem to produce very similar results. 
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The positions of the shocks are almost identical to those computed in [T9] . 
There are small deviations from their solutions in the central region along the 
axis. In this domain the flow becomes stationary with low mass density. This 



discrepancy becomes more obvious in the pressure profiles shown in Fig. |3.4 
The solid lines in these diagrams correspond to a solution obtained on a finer 
grid with a resolution of 800 x 1200 cells. In the limit r — > the pressure on 
the coarse grid is a little above 0.95, but the solution converges in both cases 
to the value stated in [T9] as the grid is refined. 

Although there is almost no difference between the two methods regarding 
the numerical results, we found that the scheme utilising the midpoint rule 
is more robust in the high Mach number regime. In simulations with Mach 
numbers above 100 we observe that the trapezoidal scheme tends to steepen 
pressure gradients which causes negative pressure in these remarkably high 
supersonic flows. 



3. 3 Detailed conservation of angular momentum 



The test introduced in this section makes use of an important property of 
inviscid rotationally symmetric flows: The coupling of angular momentum 
transport to mass transport. If one defines the specific angular momentum by 



i = Va (3.6) 



the fourth component of the system (1.7) with (3.1), (3.3) may be rewritten 
in the form 

d t [ 6 i) + V^g£ Vi ) + V n [ e £v^) = 0. (3.7) 

This equation describes the advection of angular momentum density g£ in 
curvilinear orthogonal coordinates with rotational symmetry. It is of the very 
same form as the equation for transport of mass density g, i. e. the continuity 
equation 



d t g + V^gv^ + V^QVff) = 0. (3 .1 

The time evolution of specific angular momentum is hence given by 



v* „ „ v n 



d t e = --±d i e--^d r] e. (3.9) 



tip tim 
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3. 3. 1 The mass spectrum as a constant of motion 

Following the definition of Norman et al. [20] we define the mass spectrum of 
the specific angular momentum by 



M{£) = / dm(£') 
o 



(3.10) 



with m{(!) being the mass distribution function with respect to specific angular 
momentum. For inviscid axisymmetric flows this spectrum is a constant of 
motion within a spatial region D if there is no flux across its boundary 3D 



dM 



dt 



if 



h- v\ dD = 0. 



D 



(3.11) 



n is the surface normal of the boundary and v the velocity of the fluid. With 



the help of the step function we may transform the mass integral in (3.10) 
into a volume integral 



dM 



dt 



D 



dt 



e(e-e'(Z, v ,t)) g (z, v ,t)dv(t, v ). 



D 



The spatial region D does not depend on time, hence it can be exchanged 
with the time derivative within the integration. 



dM 



dt 



D 



I d t Q{£ - £') g dV + f Q{£ - £') d t g dV (3.12) 



D 



D 



In doing so we implied that the derivative of the step function is well defined 



(see [21]). Equations (3.8) and (3.9) allow us to replace the time derivatives 
with spatial ones. For convenience we define the two-dimensional restrictions 
of the usual curvilinear differential operators by 



V ■ v = V^ + V 



r[Vf) 



and v ■ V£ = ^ dd + p- dj 

he h„ 



Therefore (3.12) becomes 



dM 



dt 



= - / qv ■ W ^r6(£ - £') dV - [ e(£ - £') V • (gv) dV. 
d J d £ J 

D D 
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The second integral may be evaluated by using integration by parts and ap- 
plying the Gaussian divergence theorem 



dt 



dM r — r 

= - gvVQ{£-£') dV - oQ{£-£')vhdA 

D dD 

+ I qv VQ(£ - £') dV = 0. 



D 



The surface integral yields zero due to the vanishing normal velocity across 
the boundary. 

The mass spectrum M(£) - although a global quantity - carries information 
about the local redistribution of angular momentum. Any kind of diffusive 
transport of the angular momentum causes a deviation from the initial spec- 
trum. Since numerical diffusion exists in any numerical scheme, we do not 
expect an exact conservation of the mass spectrum. 



3.3.2 Analytical expressions for the mass spectrum 



In general the evaluation of the integral in (3.10) is difficult and must be done 
numerically. However, for a quantitative analysis it is very useful for com- 
parison to start numerical simulations with angular momentum distributions 
for which analytical expressions of the mass spectrum exist. Since cylindrical 
coordinates {z, r, if} exhibit the natural system for the formulation of rota- 
tionally symmetric problems, we will perform all calculations in this system 
and assume symmetry with respect to tp. Let us define the surface density by 



g(r,z)dz (3.13) 



and demand that this function is well defined for r e IR + . If we claim that the 



specific angular momentum depends only on r, then the mass spectrum (3.10) 
may be written in terms of the surface density 



M(£) = 2tt / £(r)9((£ - £'{r)) r dr. (3.14) 

o 

Let us furthermore assume that £'(r) is of the form 



r 



a 



P'(r)=£ 1 — 1 with a,£ ,r eR + (3.15) 
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then the inverse function r(£') = r$ {£' /£q) is well defined and the mass 
spectrum becomes 



M{£) = 2n f T,(r)r dr. 



(3.16) 



The last integral may be evaluated analytically in some cases. We will consider 
mass distributions for which S(r) is a very simple function of r: 



(3.17) 



(1) homogeneous cylinder centred on the axis with radius R and mass M 

(2) homogeneous sphere centred on the axis with radius R and mass Mq 

2 / ,, \ 2/a\ 3/2" 



M{£) = M 



1 - 1 -Sjlt 



(3.18) 



(3) Gaussian density distribution centred on the axis with standard devia- 
tion a and mass M 



M{£) = M 



1 — exp 



I fro 

2\ a 



2/a\ 



L 



(3.19) 



These expressions become even simpler for rigid motion. In this case a = 2 
and the specific angular momentum may be expressed in terms of the constant 
angular velocity fl 



£{r) 



to 2 



r z = tt r\ 



The spectrum of a rigidly rotating sphere (|3.18|) then reduces to the formula 
given in 



3.3.3 The rigidly rotating Gaussian density distribution 

In this test we examine the disruption of a Gaussian density pulse by cen- 
trifugal forces. A dense pulse is located on the axis of symmetry within a low 
density environment. The rotational velocity in the whole computational do- 
main is initialised with constant angular velocity, i e. v v = fior 2 with Q = 10. 
All other velocity components are set to zero and the pressure is unity. The 
peak density of the pulse is £> max = 10 and the uniform density of the ambi- 
ent medium is £ m i n = 10~ 2 . The pulse is centred at (r,z) = (0,0.4) with a 
full width half mean of 0.1. On the cylindrical mesh the computational do- 
main covers a region of [0, 1] x [0, 1] with reflecting boundary conditions at all 
boundaries. In addition we switched the sign of the rotational velocity within 
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t = 0.1 



£ = 0.2 





Fig. 3.5. Schlieren type images showing the absolute value of the density gradient 
at different times; logarithmic scale between 10~ 5 (white) and 10 4 (black). The 
solution was obtained on a 800 x 800 cylindrical grid. 



the ghost zones along the axis of symmetry. All computations in this chapter 
were performed using the midpoint quadrature scheme with the monotonized 
central limiter and a Courant number of 0.4. The limiters parameter was set 
to 9 = 1.3 unless stated otherwise. 



This problem setup ensures that a huge amount of mass with low specific 
angular momentum spreads out into the whole computational domain. The 
dynamic behaviour is, at least in the beginning, mostly driven by centrifugal 
forces. Hence the dynamic time scale is dominated by the time of circulation, 
which is roughly l/fi = 0.1. Fig. 3.5 depict the time evolution of the flow up 
to t = 0.4 for a resolution of 800 x 800 cells. The forces acting on the density 
pulse accelerate the gas to supersonic speed in the radial direction forming 
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Fig. 3.6. Profile of mass density (solid Fig. 3.7. Profile of radial velocity v r 

line, left scale) and pressure (dashed (solid line) and rotational velocity v v 

line, right scale) at z = 0.4 and time (dashed line) at z = 0.4 and time 

* = 0.1. * = 0.1. 
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Fig. 3.8. Mass density and specific an- 
gular momentum (contours) at time 
t = 0.4. The scale for the contour lev- 
els is logarithmic with basis 2 starting 
at 2~ 9 . 
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Fig. 3.9. Profile of mass density (solid 
line, left scale) and specific angular mo- 
mentum (dashed line, right scale) at 
z = 0.4 and time t = 0.4. 



a bow shock. One clearly sees how the material of the dense region is driven 
outwards. Within this central region behind the shock the peak density drops 
to 5.0 at time t — 0.1 and the pressure arrives at a constant value of 0.3787 (cf. 



Fig. 3.6). Since the pressure gradient inside this bubble confined by the shock 



is zero, there are no forces acting in the vertical direction and v z remains zero. 
Furthermore the remarkable identity v r = v v holds as can be seen for z = 0.4 
in Fig. 3/7 At t = 0.2 the density pulse has flattend to a disk like structure 
causing the low pressure region to collapse roughly around t = 0.3. 
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Fig. 3.10. Mass spectrum of specific an- 
gular momentum for the initial condi- 
tion and at time t = 0.4; the resolution 
of the simulation is 800 x 800. 



Fig. 3.11. Mass spectrum of specific an- 
gular momentum at time t = 1.0 for dif- 
ferent resolutions. 



Another interesting feature forms around (r, z) = (0.5, 0.4) which is very much 
like the well known Rayleigh- Taylor instability [22.23J. The dense material of 
the pulse penetrates into the rarefied gas of the ambient medium, in this 
case driven by centrifugal forces. These instabilities grow rapidly and the flow 
becomes more and more turbulent from the time t = 0.4 (q. v. Fig. 3.13). 



The outwards driven material carries a vast amount of the mass initially con- 
centrated around (r, z) = (0,0.4). Connected to this mass flux gas with low 
specific angular momentum is transported to larger radii as can be seen in 
Fig. 3.8| and 3J3 The question under consideration is whether this redistri- 
bution in space is solely advective or partly diffusive. Therefore we compute 
the mass spectrum as described in the previous sections. In this test problem 
the exact solution can be obtained via summation of the two solutions for the 



Gaussian density pulse (Eq. 3.19) and the homogeneous cylinder (Eq. 3.17) 



with different total masses Mq. In Fig. |3. 10] this analytical solution is depicted 
in conjunction with the initial spectrum and the numerical result at time 
t = 0.4. The steep slope in this logarithmic diagram starting at £ = 10~ 6 in- 
dicates the Gaussian pulse whereas the kink at £ = 10 _1 marks the transition 
to the ambient medium. Besides some spurious data points at the lower end 
of the spectrum, the numerical result is in good agreement with the exact 
solution. Due to the limited resolution of the numerical computation, there 
are some data points with different angular momentum attached to the same 
mass. This shows up in both the initial setup and in the results for t = 0.4. 
Also there is clearly a lower limit for mass and specific angular momentum de- 
termind by the size of the grid cells and the distance between the barycentres 
of the innermost cells and the axis. 

The picture becomes completely different, if we examine the data at later 
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t = 0.6 




Fig. 3.12. Mass spectrum of specific angular momentum. Comparison of the results 
for different settings of the monotonized-central limiter at different times t; the 
resolution of the simulations is 400 x 400. 

times. After the generation of the first instabilities diffusive processes effect 
the redistribution of angular momentum. This shows up in the loss of cells 



with low mass causing a steeper slope in the spectrum seen in Fig. 3.11 Even 
more noteworthy there is a dependence on the resolution, which seems to 
be contradictory. The simulations with low resolution produce better results, 
whereas the numerical diffusion should decrease with increasing resolution. To 
study the influence of the numerical diffusion we changed the parameter of the 



monotonized-central limiter (see Eq. 3.4). In general one observes that higher 



values of 9 result in less diffusive schemes. In Fig. 3.12 the mass spectrum 



before and after the generation of the first instabilities is shown for different 
settings of the limiter. Up to t = 0.4 there is almost no dependence on the lim- 
iters parameter 9. At later times diffusive angular momentum transport occurs 
as in the previous experiment. However a comparison of the mass spectra is 
questionable if instabilities are generated. In such case a proposition regarding 
the conservation of angular momentum is of limited significance, because the 
spatial redistribution of angular momentum proceeds differently. 

To shed light on this phenomenon, we compared the vorticity of the numerical 
results at the same time for different resolutions. In cylindrical coordinates 
with rotational symmetry the components of the vorticity are given by 



(3.20) 



w z = Wr = i rv \ and w^ = 
oz r or v ' 


dv r 

dz 


dv z 
Or 


ith help of (3.9) one easily verifies, that 






d t £ = r(v z w r — v r w z ) 







(3.21) 

holds. Hence there is a close connection between the time evolution of specific 
angular momentum and the projection of the vorticity onto the r-z-plane. 
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Fig. 3.13. Absolute value of the projection of the vorticity onto the r-z- plane. Com- 
parison of the results for two different resolutions; left image: 200 x 200, right image 
800 x 800; both at time t = 1.0. 



£ = 0.4 
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Fig. 3.14. Mass spectrum of specific angular momentum. Comparison with the 
ZEUS2D code at different times t for different resolutions. 



Fig. |3.13| shows the absolute value of the projected vorticity for simulations 
performed with different resolutions. There is clearly much more turbulent 
motion in the high resolution simulation. Lots of eddies have formed and 
the turbulence seems to be much more evolved than in the results obtained 
for lower resolution. Dark regions in this diagram are an indicator for large 
vorticity and mark the surface where angular momentum is exchanged mostly 
between grid cells. A comparison by eye shows that this effective surface for 
angular momentum transport is much bigger for high resolution simulations. 
Therefore diffusive transport of angular momentum due to numerical diffusion 
becomes much more efficient. 



Finally we compare our numerical scheme to the well known second-order- 
accurate van Leer method [I] used in the ZEUS-2D code [21] developed by 
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t = 0.2 



£ = 0.4 




Fig. 3.15. Schlieren type images showing the absolute value of the density gradi- 
ent at two different times. The the numerical solution was computed on an oblate 
spheroidal mesh with a resolution of 800 x 800 grid cells; scales are identical to those 
in Fig. Boa 



Stone and Norman [25]. The ZEUS-2D code implements the numerical scheme 
for consistent angular momentum transport proposed by Norman et al. [20] 



and improved by Norman and Winkler [26]. The results depicted in Fig. 3.14 



were computed with a von Neumann and Richtmyer type scalar artificial vis- 
cosity [27]. The parameter which controls the strength was set to 2.0. For 
comparison we performed the same tests using the tensorial artificial viscos- 
ity [2H], but the results were almost identical. The mass spectra obtained for 
these simulations show diffusive transport in the low-mass and low-angular- 
momentum regime. We found that the ZEUS-2D code generates grid zones 
with zero angular momentum indicated by the horizontal branch in the left 



diagram of Fig. |3.14[ This leads to an overestimation of the mass confined in 
grid zones with low specific angular momentum. Henceforth the shape of the 
mass spectrum changes rapidly, when the first instabilities appear, especially 
for higher resolution. By contrast, the mass spectrum computed for the new 
curvilinear central-upwind scheme does not show major deviations from the 
exact solution for t < 0.6. 



3.4 The Gaussian pulse on an oblate spheroidal mesh 



The following test is meant as a demonstration for the validity of the numerical 
scheme applied to other orthogonal and rotationally symmetric coordinate 
systems. Hence we selected oblate spheroidal coordinates {£, 77, 0} (see e. g. 
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i = 0.4 

i — | — , — , — | — 

cylindrical 
oblate spheroidal 
exact 




t = 1.0 

-, 1 1 1 1 — 

cylindrical 
oblate spheroidal 
exact 



logio e lo giO l 

Fig. 3.16. Mass spectrum of specific angular momentum. Comparison of the results 
for cylindrical and oblate spheroidal grids at different times t; the resolution of the 
simulations is 800 x 800. 



chap. 21.3) with 



( x \ 



y 



( 



a^fj cos « 
a^fj sin < 



and metric scale factors 



ti£ = h v = a J (sinhc;) 2 + (sin 77)' 



a cosh £ cos 77 cos c^ 

a cosh £ cos 77 sin <; 

. a sinh £ sin rj . 



a cosh £ cos rj. 



(3.22) 



(3.23) 



The factor a is a scaling constant which is set to unity in our calculations. 
If we assume symmetry with respect to the coordinate <ft, the system of two- 



dimensional Euler equations is again given by (1.7) with flux vectors (3.1) and 



source terms (3.3). One easily proves by derivation of the scale factors that 



none of the commutator coefficients vanishes. Therefore all geometrical source 
terms have to be taken into account. 

We examine the solution for the same initial condition as described in the 
previous chapter. The Gaussian pulse with peak density p max = 10 is located 
on the axis and embedded in an ambient medium with constant density £ m ; n = 
10~ 2 . The whole computational domain is initialised with constant angular 
velocity Q = 10 and constant pressure p = 1. All boundaries are treated 
as reflecting walls. The only difference to the test setup on the cylindrical 
mesh is the shape of the computational domain with spatial extension (£, 77) G 
[0, 0.88] x [71-/8, tt/2]. As for the cylindrical mesh the resolution is 800 x 800. 



For comparison with the simulations on cylindrical grids (see Fig. 3.5) the 
numerical results for the density gradient are depicted in Fig. |3. 15 At time t = 
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0.2 the solutions on the oblate spheroidal mesh are almost identical to those 
on the cylindrical mesh. There are small differences regarding the reflection of 
shock waves at the outer boundaries. The influence of these disturbances on 
the main features becomes more important at later times (see right diagram). 
Waves reflected backwards interact with the instabilities mentioned in the 
previous chapter in a different way. Thus the solutions diverge more and more. 



This does not affect the mass spectrum in a considerable manner (see Fig 3.16 ). 
As for the cylindrical grid we observe that the spectrum is preserved very well 
up to t = 0.4. Diffusive transport at later times seems to alter the results in a 
similar way in either of the two curvilinear schemes. 



4 Conclusions 



The main advantage of central-upwind schemes is their simplicity. Apart from 
the propagation speeds of the non-linear waves no other information about 
the conservation laws under investigation is required. Therefore these schemes 
are applicable to a variety of physical problems. Although the solution of the 
Euler equations for gas dynamics is a major objective, one may also consider to 
solve the equations for ideal magneto-hydrodynamics or the Hamilton- Jacobi 
equations, (see [3"0"|31|lll] ). In this paper we extend this generality to non- 
cartesian systems, thus providing the ability to solve systems of hyperbolic 
conservation laws on orthogonal curvilinear grids. 

A computer program written in Fortran 95 has been developed in order to 
test the new numerical schemes [32] ■ So far the program is capable of solving 
the equations of inviscid gas dynamics in cartesian, polar, cylindrical, spherical 
and oblate spheroidal coordinates. In case of rotational symmetry the trans- 
port of angular momentum is included. Both quadrature rules for computation 
of fluxes and source terms are implemented. We performed excessive tests to 
verify the correctness of the numerical results. Some of them are presented in 
this paper. The solutions of two-dimensional Riemann problems discussed in 



Sections 3.1 and 3.2 are in good agreement with the cartesian results presented 



by other authors [T|19Jll6j . 

So far tests which check the detailed conservation of angular momentum are 
very rare in the literature. Norman et al. [20] perform a test similar to our 
setup, but they use a moving grid and include selfgravity in their calculations. 
The significance of these results is limited due to the fact that they study 
the mass spectrum at less than a time of circulation. Therefore we proposed 
a pure hydrodynamical test with angular momentum transport in Sec. 3.3.3| 



A comparison with the ZEUS-2D code [21] demonstrates that our numerical 
scheme leads to better results with less diffusive angular momentum transport. 
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Appendix Commutator coefficients 



For orthogonal coordinate systems all off-diagonal elements in the metric ten- 
sor vanish and the infinitesimal squared distance ds 2 is written in terms of 
metric coefficients and infinitesimal coordinate displacements according to 

ds 2 = gij dx* dx j = (ht d£) 2 + (h v dr]) 2 + (h^ dtf . (Ai) 

Hence one might define a set of orthonormal basis vectors by 

ei ~htd£ 6f, ~ h v d v e *~ h^dct> { ' 

and the corresponding set of dual 1-forms by 

& = h<: d£ £* = h v dr] ^ = h^ d0. (A3) 

Furthermore one expands all vectors, tensors and forms with respect to the 
new basis. Written in terms of orthonormal 1-forms the infinitesimal squared 
distance becomes 

ds 2 = (u£) 2 + (u^) 2 + (&*) 2 = 5 M A' 1 (A4) 

i. e. the metric coefficients are independent of the coordinates. Therefore all 
derivatives of the metric with respect to coordinates vanish and the affine 
connection is determined by 



■pfc^ 



l^icss + cm-^i) ( A5 ) 



where g kl is the inverse metric and cm denote the commutator coefficients 
defined by means of the Lie-bracket 



e h e ] 



cv? h = 9 kl c Vii gj, (A6) 



-ij ^k ~ y ^ijl °fc 



To obtain the commutator coefficients for the orthonormal basis mentioned 
above we have to apply the Lie-bracket of two basis vectors on an arbitrary 
function and carry out the derivatives. 



6p, 6fi 



_ ( 1 dh v \ 1 df ( 1 dh{\ 1 df 



v h(h v d£ J h v dr] yh^h^ dr\ ) h^ d£ 
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The non-vanishing commutator coefficients are therefore given by 







1 dh n 






1 dhcj, 


f)iv 


ivv 


h n h^ d£ 
1 dht 


c 4>i4> ~ 


c iU ~ 


1 (9/i0 


c ivi = 


~ c vii = 


h^h v di] 

1 dfl£ 


4>i)4> ~~ 


f/4>4> ~ 


h^hr, drj 
1 dh v 


c iU = 


~ c 4>ii = 


h^h^ (90 


fjcjyf) ~~ 


4>fjv ~ 


hjji^ d(j) 



(A7) 
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